
(* 
**********************************************************

		 COSY_PAK package chap4.m		

**********************************************************
*)


If [ TrueQ[ $VersionNumber >= 2.0 ],
     Off[ General::spell ];
     Off[ General::spell1 ] ]
     
(*  B E G I N     P A C K A G E  *)


BeginPackage["COSYPAK`chap4`"]	     

Routh::usage = " ";


Begin["`Private`"]

Routh[charpoly_,s_,z_:zero] :=
Module[{Clist,ncoeff,list1,list2,Routhtable,nnn,temp,jj,i,first,k},
	Clist = Reverse[CoefficientList[charpoly,s]];
	ncoeff = Length[Clist];
	nnn = ncoeff;
	
	If[OddQ[ncoeff], ncoeff = ncoeff +1;   (* make Coefflist to be *) 
	AppendTo[Clist,0]];  			(* even number	*) 
	list1 = {}; list2= {}; Routhtable = {};
	Do[ AppendTo[list1,Clist[[2 k -1]]],{k,1,ncoeff/2}];	
	Do[ AppendTo[list2,Clist[[2 k]]],{k,1,ncoeff/2}];
	If[list2[[1]] == 0 , list2[[1]] = z ];
	
	Routhtable = Join[{list1},{list2}];

	
	For[jj = nnn -2, jj > 0 , jj--,
	    temp={};
		Do[elm=Simplify[(list2[[1]] list1[[i+1]]-list1[[1]] list2[[i+1]])/
						list2[[1]] ];
		   If[ (i == 1) && Abs[elm] < 0.00001 , elm= z ];		
		   AppendTo[temp,elm],
		  {i,1,ncoeff/2 - 1}];
		AppendTo[temp,0];	  		
		Routhtable = Append[Routhtable,temp];
		list1 = list2; list2 = temp
		];
		
	first = Table[Routhtable[[i,1]],{i,1,Length[Routhtable]}];	
	Print["The Routh's Table:"];
	Print[" "];
	Print[MatrixForm[Routhtable]];
	Print[" "];
	Print["The first column of the Routh's table is: "];
	Print[" "];
	Print[ first ];
	Print[" "];
	Do[ If[first[[i]] < 0, Print["Sign Change!! Unstable root(s) exist."]],
		{i,1,Length[first]}];
	Return[Routhtable]	 	
	]/;PolynomialQ[charpoly,s]

End[];
EndPackage[];

(*  E N D     P A C K A G E  *)

If [ TrueQ[ $VersionNumber >= 2.0 ],
     On[ General::spell ];
     On[ General::spell1 ] ]

